Modelling the Impact of NETosis During the Initial Stage of Systemic Lupus Erythematosus

The development of autoimmune diseases often takes years before clinical symptoms become detectable. We propose a mathematical model for the immune response during the initial stage of Systemic Lupus Erythematosus which models the process of aberrant apoptosis and activation of macrophages and neutrophils. NETosis is a type of cell death characterised by the release of neutrophil extracellular traps, or NETs, containing material from the neutrophil’s nucleus, in response to a pathogenic stimulus. This process is hypothesised to contribute to the development of autoimmunogenicity in SLE. The aim of this work is to study how NETosis contributes to the establishment of persistent autoantigen production by analysing the steady states and the asymptotic dynamics of the model by numerical experiment. Supplementary Information The online version contains supplementary material available at 10.1007/s11538-024-01291-3.

Parameter values used in the numerical experiments in the main text: 2 Estimates for the analytical solution The estimates for non-negativity and boundedness of the solutions of the model are made for times t > 0. The model equations are Proposition 1 The solutions to (1) remain non-negative in time t > 0 for non-negative initial values x 1 (0), x 2 (0), y(0), z(0).
Proof.Since the right-hand-side of ( 1) is continuous and continuously differentiable, the initial value problem has a unique solution for some positive time.We have that Assume that the solution of (1) eventually leaves the non-negative orthant of R 4 and take τ := inf{t > 0 : Consider equation (1c) and observe that for t ≤ τ , y ′ ≥ −(µ 4 + µ 5 y)y.

Parameter estimation
Antigen: The value range for the yield of autoantigen from NETosis α is chosen between 0 and 25% of a single cell mass, or up to 0.25 ng/cell, if the mass of a neutrophil cell is taken as 1 ng (10 −9 g), in line with the estimate of cell mass made in [2].
The value ν 1 = 0.5 corresponds to the observation from in vitro measurement of the time for leakage of nucleosome content from apoptotic Jurkat cells [3].This study shows that apoptotic cells transition to necrotic matter by release of nucleosomes in 24-48 hours (hence the rate ν 1 = 0.5/day).However, since this estimate comes from an in vitro assay, we also take a 10-fold lower value (ν 1 = 0.05/day) and vary ν 1 as a bifurcation parameter to account for the action of any additional innate immune mechanisms that we do not model explicitly.In accordance with this estimate we set the removal of apoptotic matter and necrotic matter to occur on a much faster scale of 2-2 1  2 hours leading to the estimates of rates µ 1 , µ 2 ∈ [10, 12]/day.

Macrophages:
The removal rate of macrophages µ 4 = 0.2 corresponds to a life-span of 5 days, comparable to the estimate provided in [4] (life-span of 4-7 days for blood monocytes).The values for β 1 , β 2 , σ 2 , µ 5 are chosen to correspond to a maximum density of recruited macrophages in the order of 10 6 − 10 7 cells/ml, similar to the estimates of [1].
The parameter σ 3 is chosen to bring the steady state value of z = σ 3 /µ 3 for the normal state E 0 within the range of neutrophil count for healthy volunteers 3×10 6 to 10 7 cells/ml [6].We let accordingly σ 3 take values in the interval (3.7×10 6 , 5×10 6 ).

Additional parameter sets
We also perform simulations for modified rates of NETosis and neutrophil dynamics using the following parameter set: (P.9)These values meet the conditions for multistationarity of the model ( 1) for α = 0 described in the main text, and the existence of bistability in the vicinity of α = 0 is verified numerically.
Figures S.4 and S.5 show the bifurcation structure within the biologically relevant range of α.A supercritical Hopf bifurcation H + on the blue coexistence branch at α = 2.91×10 −4 leads to the emergence of a stable limit cycle.We note that for much larger values of α in the interval [0.0122, 0.0125] we observe multistationarity on the blue branch, with three coexistence-type states which exist in parallel and all of them being unstable for α ∈ [0.0122, 0.0123] (not shown).However, this interval is outside the biologically assumed range for α.
We also explore the effect of varying the rates σ x 1 , for parameter values given in (P.9) and α ∈ (0, 3.5×10 −4 ).Zoomed panels for α ∈ (0, 8×10 −6 ) and (0, 6×10 −4 ) show the bifurcation structure in more detail.The branch of states E 0 is shown in green, the branch of E 1 in magenta, the two branches of type E * in blue and orange.The black dots represent the branching points between E 0 , E 1 , E * , and the red dot -the supercritical Hopf bifurcation.The numerical bifurcation analysis presented in Figures S.12 and S.13 reveals a regime where a stable limit cycle coexists with a steady state of type E * (blue branch).Another supercritical Hopf bifurcation H ++ makes the blue coexistence branch lose stability at α ≈ 0.282×10 −4 .The oscillations inside the limit cycle gradually increase in amplitude.
If we keep the parameters (P.10) but reduce ν 2 to ν 2 = 0.033, we obtain similar dynamics, except that the transcritical bifurcations are shifted for larger values of α.The numerical bifurcation analysis presented in Figures S.14 and S.15 reveals a regime where a stable limit cycle again coexists with a steady state of type E * (blue branch).The blue coexistence-type branch E * is locally asymptotically stable over the entire biologically relevant range of α before a supercritical Hopf bifurcation occurs at α = 2.8195×10

1 (Figure S. 4 :
Figure S.4: Bifurcation diagram α vs.x 1 , for parameter values given in (P.9) and α ∈ (0, 3.5×10 −4 ).Zoomed panels for α ∈ (0, 8×10 −6 ) and (0, 6×10 −4 ) show the bifurcation structure in more detail.The branch of states E 0 is shown in green, the branch of E 1 in magenta, the two branches of type E * in blue and orange.The black dots represent the branching points between E 0 , E 1 , E * , and the red dot -the supercritical Hopf bifurcation.

Figure S. 12 :
Figure S.12: Bifurcation diagram, α vs. x 1 for α ∈ (0, 1.2×10 −4 ).The branch of state E 0 is plotted in green, branch of E * in blue.The red lines mark the minima and maxima values in the limit cycle which arises from the supercritical Hopf bifurcations at H − , H + (parameter values in (P.10)).
Figure S.17: Fractions of the two types of antigen as function of β 1 , β 2 in the stable steady state of type E * .Remaining parameter values given in (P.11).
Figure S.14: Bifurcation diagram, α vs. x 1 for α ∈ (0, 1.2×10 −4).The branch of state E 0 is plotted in green, branch of E * in blue.The red lines mark the minima and maxima values in the limit cycle which arises from the supercritical Hopf bifurcations at H − , H + (parameter values in (P.10), but ν 2 = 0.033).